# ============================================================================
# Large Scale GPU Performance Test
# OpenSees Tcl Script for GPU vs CPU Performance Comparison
# ============================================================================

wipe
puts "Large Scale GPU Performance Test Started"

# ============================================================================
# Test Configuration - Larger Problem Size
# ============================================================================

set earthquakeFile "GroundMotions/earthquake.txt"
set timeStep 0.001
set analysisDuration 2.0      ;# Short analysis for performance test

set timestamp [clock format [clock seconds] -format "%Y%m%d_%H%M%S"]
set outputDir "Results/PerformanceTest_$timestamp"
file mkdir $outputDir

# ============================================================================
# Large Model Creation - 4x larger than original
# ============================================================================

model basic -ndm 2 -ndf 2

# Geometry parameters - Increased size
set soilWidth    20.0    
set soilHeight   10.0    
set nEleX        200      ;# 2x increase (was 20)
set nEleY        100      ;# 2x increase (was 10)

# This gives us: 40x20 = 800 elements (vs 200 original)
# Nodes: 41x21 = 861 nodes (vs 231 original)  
# DOFs: ~1722 DOFs (vs 462 original)

set dx [expr $soilWidth/$nEleX]
set dy [expr $soilHeight/$nEleY]

# Material parameters
set rho        1800.0   
set E          50.0e6   
set nu         0.3      

puts "Large model: ${nEleX}x${nEleY} = [expr $nEleX*$nEleY] elements"
puts "DOFs: approximately [expr ($nEleX+1)*($nEleY+1)*2] DOFs"

# Create nodes
set nodeTag 1
for {set j 0} {$j <= $nEleY} {incr j} {
    for {set i 0} {$i <= $nEleX} {incr i} {
        set x [expr $i * $dx]
        set y [expr $j * $dy]
        node $nodeTag $x $y
        incr nodeTag
    }
}

# Define material
nDMaterial ElasticIsotropic 1 $E $nu $rho

# Create elements
set eleTag 1
for {set j 0} {$j < $nEleY} {incr j} {
    for {set i 0} {$i < $nEleX} {incr i} {
        set n1 [expr $j*($nEleX+1) + $i + 1]
        set n2 [expr $j*($nEleX+1) + $i + 2]  
        set n3 [expr ($j+1)*($nEleX+1) + $i + 2]
        set n4 [expr ($j+1)*($nEleX+1) + $i + 1]
        
        element quad $eleTag $n1 $n2 $n3 $n4 1.0 "PlaneStrain" 1
        incr eleTag
    }
}

# Boundary conditions
for {set i 1} {$i <= [expr $nEleX+1]} {incr i} {
    fix $i 1 1
}

for {set j 1} {$j <= $nEleY} {incr j} {
    set leftNode [expr $j*($nEleX+1) + 1]
    fix $leftNode 1 0
    set rightNode [expr $j*($nEleX+1) + ($nEleX+1)]
    fix $rightNode 1 0
}

# ============================================================================
# Quick Gravity Analysis
# ============================================================================

puts "Performing gravity analysis..."

pattern Plain 1 Linear {
    for {set j 0} {$j <= $nEleY} {incr j} {
        for {set i 0} {$i <= $nEleX} {incr i} {
            set nodeNum [expr $j*($nEleX+1) + $i + 1]
            if {$j == 0} continue
            
            set tributaryArea [expr $dx * $dy]
            if {($i == 0 || $i == $nEleX)} {set tributaryArea [expr $tributaryArea * 0.5]}
            if {($j == $nEleY)} {set tributaryArea [expr $tributaryArea * 0.5]}
            if {($i == 0 || $i == $nEleX) && ($j == $nEleY)} {set tributaryArea [expr $tributaryArea * 0.5]}
            
            set bodyForce [expr -$rho * 9.81 * $tributaryArea]
            if {$tributaryArea > 0.0} {
                load $nodeNum 0.0 $bodyForce
            }
        }
    }
}

constraints Plain
numberer Plain
system BandGeneral
test NormDispIncr 1.0e-6 25
algorithm Newton
integrator LoadControl 0.1
analysis Static

analyze 10
loadConst -time 0.0
wipeAnalysis

# ============================================================================
# Create Simple Earthquake Loading
# ============================================================================

# Create simple sinusoidal ground motion for testing
timeSeries Sine 1 0.0 $analysisDuration 1.0 -factor [expr 0.1*9.8] 
pattern UniformExcitation 2 1 -accel 1

# Setup damping
set dampRatio 0.05
set omega1 [expr 2*3.14159*1.0]   
set omega2 [expr 2*3.14159*10.0]  
set a0 [expr $dampRatio*2*$omega1*$omega2/($omega1+$omega2)]
set a1 [expr $dampRatio*2/($omega1+$omega2)]
rayleigh $a0 0.0 0.0 $a1

# ============================================================================
# Performance Comparison Test
# ============================================================================

set testSteps 50  ;# Number of steps for performance test

puts "=== Starting Performance Test ==="
puts "Test duration: [expr $testSteps * $timeStep] seconds ($testSteps steps)"

# Test BandGeneral
puts "Testing BandGeneral solver..."
constraints Plain
numberer Plain
system BandGeneral
test NormDispIncr 1.0e-6 25
algorithm Newton
integrator Newmark 0.5 0.25
analysis Transient

set bandStartTime [clock clicks -milliseconds]
set ok 0
for {set step 1} {$step <= $testSteps} {incr step} {
    set ok [analyze 1 $timeStep]
    if {$ok != 0} {
        puts "BandGeneral failed at step $step"
        break
    }
}
set bandEndTime [clock clicks -milliseconds]
set bandTime [expr $bandEndTime - $bandStartTime]
puts "BandGeneral: $testSteps steps completed in ${bandTime} milliseconds"

# Clear and test GPUSolver
wipeAnalysis

puts "Testing GPUSolver..."
constraints Plain
numberer Plain
system GPUSolver
test NormDispIncr 1.0e-6 25
algorithm Newton
integrator Newmark 0.5 0.25
analysis Transient

set gpuStartTime [clock clicks -milliseconds]
set ok 0
for {set step 1} {$step <= $testSteps} {incr step} {
    set ok [analyze 1 $timeStep]
    if {$ok != 0} {
        puts "GPUSolver failed at step $step"
        break
    }
}
set gpuEndTime [clock clicks -milliseconds]
set gpuTime [expr $gpuEndTime - $gpuStartTime]
puts "GPUSolver: $testSteps steps completed in ${gpuTime} milliseconds"

# ============================================================================
# Results Analysis
# ============================================================================

puts "\n=== PERFORMANCE COMPARISON RESULTS ==="
puts "Model size: [expr $nEleX*$nEleY] elements, [expr ($nEleX+1)*($nEleY+1)*2] DOFs"
puts "Test steps: $testSteps"
puts ""
puts "BandGeneral time: ${bandTime} ms"
puts "GPUSolver time:   ${gpuTime} ms"

if {$bandTime > 0 && $gpuTime > 0} {
    set speedup [expr double($bandTime) / $gpuTime]
    puts ""
    if {$speedup > 1.0} {
        puts "GPU SPEEDUP: [format %.2f $speedup]x faster than CPU"
        puts "GPU time savings: [expr $bandTime - $gpuTime] ms ([format %.1f [expr ($bandTime-$gpuTime)*100.0/$bandTime]]%)"
    } else {
        set slowdown [expr double($gpuTime) / $bandTime]
        puts "GPU SLOWDOWN: [format %.2f $slowdown]x slower than CPU"
        puts "Possible reasons for poor GPU performance:"
        puts "- Problem size too small for GPU architecture"
        puts "- GPU memory transfer overhead"
        puts "- GPU implementation not fully optimized"
        puts "- Single-precision vs double-precision differences"
    }
}

# Save results
set resultFile [open "$outputDir/performance_results.txt" w]
puts $resultFile "OpenSees GPU Performance Test Results"
puts $resultFile "Model: ${nEleX}x${nEleY} elements, [expr ($nEleX+1)*($nEleY+1)*2] DOFs"
puts $resultFile "Test steps: $testSteps"
puts $resultFile ""
puts $resultFile "BandGeneral time: ${bandTime} ms"
puts $resultFile "GPUSolver time: ${gpuTime} ms"
if {$bandTime > 0 && $gpuTime > 0} {
    set speedup [expr double($bandTime) / $gpuTime]
    puts $resultFile "Speedup: [format %.3f $speedup]x"
}
close $resultFile

puts "\nResults saved to: $outputDir/performance_results.txt"
puts "Performance test completed."